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Identification  of  Partially  Obscured  Objects  in  Two  and  Three 
Dimensions  by  Matching  of  Noisy  'Characteristic  Curves' 

Jacob  T.  Schwartz  and  Micha  Sharir 

Computer  Science  Department 

Courant  Institute  of  Mathematical  Sciences 

New  York  University 

ABSTRACT 

The  problem  of  matching  partially  obscured  noise-corrupted  images  of 
composite  scenes  in  two  and  three  dimensions  is  analyzed.  We  describe 
efficient  methods  for  smoothing  the  noisy  data  and  for  matching  portions  of 
the  observed  object  boundaries  (or  of  characteristic  curves  lying  on  bounding 
surfaces  of  3-D  objects)  to  pre-stored  models.  Initial  experiments  showing  the 
efficacy  of  this  procedure  are  reported. 

1.   Introduction 

The  'object  identification'  or  'model-based  vision'  problem  can  be  defined  as  follows: 
We  are  given  an  image  of  an  object  0,  which  is  assumed  to  be  a  member  of  a  collection  of 

possible  objects  0Y On  known  a  priori  (e.g.  Ox,  .  .  .  ,On  may  be  a  collection  of  parts 

involved  in  an  assembly/disassembly   process,   in  an  industrial  environment  in  which  only 

objects  drawn  from  the  set  Ox 0„  can  legitimately  appear.)    However,  the  object  O  is 

observed  from  an  unknown  angle,  and  moreover  only  part  of  O  is  observed,  the  rest  of  its 
surface  being  'obscured'  in  some  manner  (usually  by  other  objects  belonging  to  the  family 

Oj On  which  happen  to  overlap  0;  see  below).  Our  problem  is  to  identify  the  object  O 

as  one  of  the  elements  in  the  list  of  'model'  objects  Ox 0„  given  a  priori  and  to 


determine  the  position  and  orientation  at  which  O  is  being  viewed. 

This  is  the  object  identification  problem  for  a  single  object.  We  may  also  wish  to 
consider  the  more  difficult  identification  problem  for  compound  scenes  ('bin  image' 
problem),  in  which  the  image  to  be  analyzed  is  allowed  to  contain  several  objects 
0,  O' ,0"  •  •  •  all  drawn  from  the  specified  collection  Oit  .  .  .  ,0„  of  model  objects,  but 
viewed  from  an  angle  in  which  the  images  of  O,  O' ,  O"  ■  •  ■  can  overlap  and  can  obscure 
parts  of  each  other. 

This  problem  is  fundamental  to  computer  vision  and  has  been  discussed  repeatedly  in 
the  literature;  see  [BB82]  for  a  review  of  standard  approaches.  The  approach  to  be  described 
differs  from  the  feature-based  techniques  commonly  used  (see,  e.g.  [BC82]  and  [BHH84])  in 
that  it  matches  observed  objects  directly  to  models  but  in  a  more  robust  and  efficient  manner 
than  other  related  schemes  that  have  been  reported  (see  [ABBF84],  [KK85],  [TMV85])  in 
that  it  makes  only  very  limited  use  of  tangent  direction  estimates  and  uses  high-speed 
algorithms  to  calculate  the  rotation  and  translation  which  best  match  an  observed  object  to  a 
prestored  model.  The  approach  used  bears  some  relation  to  that  described  in  [GLP84], 
[GLP85].  but  is  independent  of  any  assumption  that  the  objects  to  be  identified  should  be 
polygonal  or  polyhedral. 

The  present  paper  will  consider  the  identification  problem  largely  in  two  dimensions, 
but  with  some  initial  attention  to  the  somewhat  more  difficult  three-dimensional  case.  In 
both  two  and  three  dimensions  we  will  make  the  simplifying  assumption  that  images  which 
represent  the  observed  objects  0,  0' ,  0"  ■  ■  •  in  their  true  geometric  size  are  available.  In 
the  2-D  case,  this  would  be  the  case  for  objects  cut  out  of  cardboard  or  sheet  metal,  observed 
by  a  camera  a  known  fixed  distance  above  them.  In  the  3-D  case,  the  'fixed-scale'  images 
assumed  can  be  produced  with  the  aid  of  a  'depth'  sensor  which  produces  an  'image'  showing 
the  true  3-D  position  of  each  point  in  every  observed  body  surface.  Camera-based  depth 
sensors  of  this  kind  are  now  becoming  common;  see  [S83,  BS84]  for  an  account  of  one  such 
sensor. 


Through  Section  6  below  discusses  surface-based  methods  for  object  identification;  the 
principal  approach  to  the  object  identification  problem  that  will  be  outlined  in  the  present 
paper  is  based  upon  a  notion  of  'rigidly  embedded  curve'.  This  is  any  curve  C  defined  by  the 
geometry  of  an  object  O  and  fixed  in  0,  so  that  when  O  undergoes  any  rotation  and 
translation  (i.e.  Euclidean  motion)  the  curve  C  moves  with  it.  For  2-D  objects,  these  curves 
can  simply  be  the  object  boundaries;  for  3-D  objects,  other  somewhat  more  subtly  defined 
curves,  for  example  curves  of  specified  constant  principal  curvature  on  the  observed  surface, 
or  lines  of  maximum  curvature  ('ridges')  can  be  used  instead.  (A  somewhat  more  detailed 
discussion  of  various  potentially  useful  characteristic  curves  "for  3-D  object  surfaces  is  found 
in  Section  4  below.) 

The  main  property  that  we  require  these  curves  to  possess  is  that  they  should  be 
relatively  immune  to  noise.  That  is,  if  we  take  a  model  object  O  with  a  characteristic  curve  C 
on  it,  and  then  perturb  O  by  seme  random  noise  of  maximum  size  t  to  obtain  a  'noisy'  image 
O'  of  it,  then  the  corresponding  curve  C  on  0'  should  not  deviate  from  C  by  more  than 
some  function  of  e  which  tends  to  zero  with  c.  As  we  shall  see  below,  to  enforce  such 
behavior,  some  preliminary  'smoothing'  will  be  needed  before  C"  is  calculated.  More  details 
concerning  a  smoothing  technique  which  seems  appropriate  are  given  in  Section  2  below. 

Once  we  have  defined  the  curves  to  be  used,  our  proposed  identification  method 
proceeds  as  follows: 

(a)  The  characteristic  curves  C  on  the  observed  object  and  on  all  the  model  objects  are 
parametrized  by  arc  length  s,  giving  (vector  valued)  functions  c(s)  (for  the  observed  object) 
and  ct(s)  (for  the  various  model  objects). 

(b)  We  find  the  offset  s0,  the  index  /',  and  the  Euclidean  motion  E  which  cause  the  function 
Ec(s)  to  match  the  model  curve  c,(j  +  j0)  most  closely  (the  offset  s0  is  needed  in  cases  in 
which,  because  of  obscuration,  only  a  portion  of  the  curve  C  can  be  obtained  from  the 
viewed  image).    The  orientation  of  O'  (relative  to  the  model  object  O  which  it  matches)  is 


then  defined  by  the  transformation  E.  We  shall  show  that  determination  of  E  and  s0  can  be 
accomplished  rapidly  if  we  measure  the  distance  from  c(s)  to  c,(s  +  s0)  in  a  least-squares 
sense.  Indeed,  this  mode  of  approximation  will  be  seen  to  allow  £  to  be  calculated  explicitly 
and  s0  to  be  determined  by  a  fast  Fourier-transform  based  method,  and  with  the  happy 
consequence  that  the  total  cost  of  matching  is  only  0(nm  log  m),  n  being  the  total  number  of 
model  objects  (or  more  precisely  the  total  number  of  characteristic  curves),  and  m  the 
number  of  (evenly-spaced)  points  of  interpolation  introduced  along  each  curve  to 
approximate  it.  (A  hashing  technique  for  removing  the  factor  n  from  this  last  cost  formula  is 
described  in  the  companion  paper  [KSSS85].) 

(c)  As  already  noted,  observation  of  any  characteristic  curve  C  of  the  viewed  image  O'  of 
some  model  object  O  will  always  produce  a  variant  C*  of  C  that  is  perturbed  by  addition  of  a 
certain  amount  of  noise.  Since  this  noise  can  introduce  irregularities  in  C  that  distort  its  arc 
length  significantly,  matching  must  be  preceded  by  a  smoothing  operation  which  takes  C*  and 
brings  it  back  closer  to  the  original  curve  C. 

We  now  proceed  to  explain  how  all  these  various  steps  can  be  accomplished. 

2.   Smoothing  Noise-Distorted  Curves  by  Finding  Shortest  Local  Approximations 

Existing  methods  for  smoothing  noise-corrupted  curves  ordinarily  use  some  band-pass 
convolution  operation.  However,  we  will  argue  that  in  preparing  for  curve-matching  a  rather 
different  approach,  which  smooths  such  curves  by  finding  shortest  paths  lying  in  the 
neighborhood  of  an  observed  noisy  curve,  can  be  preferable.  To  this  end,  we  begin  by 
giving  the  following  simple 

Definition  2.1:  Let  c  be  a  parametrized  curve  in  Euclidean  n-space  E",  and  let  L  be  its 
length.   Then 

(i)  The  e-neighborhood  Ut(c)  of  c  is  the  set  of  all  points  in  E"  whose  distance  from  some 
point  of  c  is  at  most  e. 
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(ii)  The  curve  c  is  (e,8)-stable  if  every  curve  in  Ut(c)  which  connects  the  endpoints  of  c  has 
length  at  least  L  -  8. 

It  is  clear  that  every  simple  (i.e.  non  self-intersecting)  smooth  curve  is  (0,0)-stable,  and 
thus  will  remain  (e,B)-stable  for  every  positive  8  and  sufficiently  small  e. 

Suppose  next  that  the  curve  c*  is  obtained  from  c  by  addition  of  some  perturbing 
(noise)  function  of  maximum  magnitude  of  at  most  e,  so  that  c*  lies  in  Ut(c)  and  connects 
two  points  a,  b  that  lie  within  €  of  the  corresponding  endpoints  of  c.  We  shall  suppose  that 
all  the  curves  with  which  we  deal,  here  and  below,  are  parametrized  in  terms  of  their  arc 
length  s,  and  that  c*  is  differentiable.  To  smooth  c*,  what  we  want  to  do  is  to  derive 
another  parametrized  curve  c' (s)  from  it  such  that  \c(s)-c'(s)\  is  uniformly  small  over  the 
whole  length  of  c;  once  this  has  been  done,  it  becomes  possible  to  use  the  matching 
procedure  sketched  in  (b)  above  to  match  the  curve  c'  to  the  model  curve  c. 

For  this  we  define  c'  to  be  the  shortest  curve  in  Ut(c*)  connecting  the  endpoints  of  c*. 
To  show  that  c'  has  the  required  property,  we  can  reason  as  follows.  Plainly  c'  lies  within 
U2t(c)  and  c  lies  within  Ut(c*).  Suppose  that  c  is  (2t,8)-stable,  and  let  V  (resp.  L)  be  the 
length  of  c'  (resp.  c).  Then  by  definition  t'2i-8.  Moreover,  the  endpoints  a,  b  of  c* 
can  be  connected  by  a  curve  which  goes  from  a  to  its  nearest  point  on  c,  then  along  c  to  a 
point  nearest  to  b,  and  finally  from  there  to  b.  Since  c  lies  within  Ut(c*)  and  c'  is  the 
shortest  curve  in  Ut(c*)  with  the  same  endpoints,  we  must  have  L  +  2i2i'.  Thus 
\L  -  L'\s  max(8,2e). 

Take  a  point  c'(s)  along  c' .  This  plainly  lies  within  a  distance  2t  of  some  point  c(s,)  of 
c.  Note  that  we  can  connect  the  initial  point  of  c  to  its  final  point  by  going  along  c  to  cisj, 
then  along  a  straight  segment  to  c'(s),  then  along  c'  to  a  point  nearest  to  the  final  point  b  of 
c,  and  finally  along  a  straight  segment  to  b.  This  whole  curve  plainly  lies  within  U2t(c)  and 
its  length  is  at  most 

sx  +  V  -  s  +  4t  «;  L  +  jj  -  s  +  6e  . 
Thus  we  must  have  L  -  8  ■<•  L  +  6«  +(sx  -  s),  i.e.  s  -  sx  s  6t  +  8.    Similarly,  we  can 


connect  the  final  point  of  c  to  its  initial  point  by  going  back  along  c  to  c(sj,  then  along  a 
straight  edge  to  c'(s),  then  to  the  point  on  c'  nearest  to  the  initial  point  a  of  c,  and  finally 
along  a  straight  edge  to  a.  Arguing  as  before  we  find  that  i-5st  +  4t+(j-  sj,  i.e. 
i:  -  i  s  4e  +  8.  This  shows  that  \s  -  s^  rs  6t  +  8,  from  which  the  following  lemma 
follows  at  once. 

Lemma  2.2:  If  the  curve  c(s)  is  (2e,8)-stable,  c*  is  derived  from  c  by  adding  a  noise  signal 
of  maximum  modulus  at  most  e,  c'  is  the  shortest  curve  within  Ut(c')  connecting  the 
endpoints  of  c*,  and  all  curves  are  parametrized  by  their  arc-length,  then 
\c'(s)  -  c(s)\  <  8t  +  8  for  all  s  in  the  common  domain  of  c  and  c' .  Moreover  the  lengths  of 
c  and  of  c'  differ  by  at  most  max  (8,2c). 

This  lemma  shows  that,  if  c  is  stable  in  me  sense  of  Definition  2.1,  with  sufficiently 
small  e  and  8,  the  curve  c'  must  match  c  well  in  the  sense  that  we  require. 

To  be  sure  that  c'  can  be  obtained  efficiently  from  the  observed  curve  c*  we  must  now 
show  how  to  calculate  shortest  paths  lying  in  the  neighborhood  of  a  given  noisy  polygonal 
curve  c*  and  connecting  the  endpoints  of  c*. 

An  algorithm  for  finding  shortest  paths  near  a  given  polygonal  curve  In  two  and  more 
dimensions. 

We  first  consider  the  relatively  simple  case  of  plane  curves,  and  begin  by  considering 
the  restricted  case  in  which  c  is  a  polygonal  (plane)  curve  which  is  the  graph  of  a  piecewise 
linear  function  y  =  c(x),  x  €  [a,b].  Let  c^x),  c2(x)  be  two  additional  piecewise  linear 
functions  satisfying  c^x)  <  c(x)  <  c2(x)  for  each  a  <  x  s  b,  and  let  U  denote  the  closed 
polygonal  strip  between  c{  and  c;,  i.e. 

U  -  {(x,y)  :  a  =s  x  =s  b,   cx(x)  =s  y  s  c2(x)}  . 
(Note  that  in  smoothing  a  polygonal  graph  c*  as  in  the  preceding  section,  we  can  take 

c^x)  =  c*(x)  -  €,      and     c2(x)  =  c'(x)  +  t,     x  €  [a,b].)     Let     A  =  (a,c(a))     and     let 

B  =  (b,c(b))  be  the  initial  and  final  points  of  c.  Our  goal  is  to  find  the  shortest  (necessarily 
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polygonal)  path  p  within  U  connecting  A  to  B .  The  technique  presented  below  achieves  this 
in  time  linear  in  the  number  of  corners  of  U.  It  is  similar  to  a  more  general  but  more 
complex  linear-time  technique  for  finding  shortest  paths  with  simple  polygons  due  to  Lee  and 
Preparata  [LP84]. 

For  each  point  r  6  U  let  it(z)  denote  the  shortest  path  from  A  to  z  within  U.  It  is  clear 
that  in  the  restricted  situation  considered  each  of  the  paths  it(z)  must  be  a  polygonal  path 
which  is  the  graph  of  a  single- valued  function  of  x,  whose  corners  must  lie  in  the  set  N 
consisting  of  the  points  A  and  B,  together  with  all  corners  of  the  graphs  of  clt  c2.  For  each 
point  z  €  U  let  v(z)  €  N  denote  the  initial  endpoint  of  the  last  straight  segment  in  it(z).  To 
obtain  p  we  will  compute  the  map  v  for  all  corners  in  N  by  using  the  following  left-to-right 
sweeping  technique.  For  each  t  €  [a, b]  let  /(f)  denote  the  vertical  line  x  =  r,  and  let  Q(t) 
denote  a  list  of  disjoint  open  intervals  7: Ij  along  /(r)  having  the  following  properties: 

(i)  The  union  of  the  closures  of  the  intervals  in  Q(t)  contains  /(r)  D  U,  and  in  the  list  Q(t) 
the  intervals  Iit...,Ij  occur  in  ascending  vertical  order. 

(ii)  For  each  interval  Ik  there  exists  a  unique  point  w  €  jV  such  that  for  all  z  6  Ik  we  have 
v(z)  =  w.  Let  us  agree  to  designate  this  property  by  saying  that  Ik  is  contained  in  the  zone  of 
influence  of  the  corner  w. 

Let  I  be  one  of  the  intervals  in  Q(t),  and  assume  that  it  is  contained  in  the  zone  of 
influence  of  some  w  €  N.  Let  Ll  and  L2  be  the  two  lines  connecting  w  to  the  two  endpoints  of 
7.  Then,  if  t'  is  slightly  larger  than  t  and  no  corner  in  N  has  an  abscissa  lying  in  [t,t'],  it  is 
easily  checked  that  the  intersection  of  l(t')  with  the  wedge  bounded  between  Lx  and  L:  will 
yield  an  interval  /'  in  Q{t')  which  is  contained  in  the  zone  of  influence  of  w,  and  furthermore 
all  intervals  in  Q(t')  are  obtained  in  a  similar  manner  from  the  intervals  in  Q(t).  In  other 
words,  if  we  let  l(t)  move  rightward,  then,  as  long  as  it  does  not  pass  through  a  corner  in  N, 
the  combinatorial  structure  of  Q{t)  remains  unchanged,  and  the  endpoints  of  each  interval  / 
in  Q(t)  simply  move  along  straight  rays  which  emerge  from  the  point  w  e  N  in  whose  zone 
of  influence  /  lies. 
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Next  suppose  that  for  some  f0,  the  line  /(f0)  does  pass  through  a  point  z  €  N,  and 
without  loss  of  generality  suppose  that  z  belongs  to  the  graph  of  cv  Suppose  that  just  before 
reaching  f0,  the  list  Q(i)  consisted  of  the  intervals  7lf  .  .  .  ,7m.  Let  w  =  v(z),  and  let  lj  be  that 
interval  which  is  contained  in  the  zone  of  influence  of  w.  Then  if  /(f)  is  moved  slightly  to  the 
right  of  r0,  the  list  Q(t)  will  have  one  of  the  following  two  forms: 

(i)  Let  e  be  the  straight  segment  in  the  graph  of  cl  whose  left  endpoint  is  z.  If  e  lies  below 
the  line  L  passing  through  w  and  z,  then  the  list  Q  will  consist  of  the  intervals 
I0,I*,Ij+i.  .  .  .  ,/,„,  where  70  is  the  intersection  of  /(f)  with  the  triangular  wedge  (having  z  as 
apex)  bounded  by  e  and  L,  and  I*  is  the  intersection  of  /(f)  with  the  wedge  bounded  by  L 
and  the  line  connecting  w  to  the  upper  endpoint  of  /,.  In  other  words,  a  lower  portion  of  lj  is 
truncated  at  z,  a  new  interval  70  is  added  below  it,  and  the  intervals  lx,  .  .  .  ,L^l  below  lj  are 
all  removed  from  the  list  Q.  Note  that  70  belongs  to  the  (new)  zone  of  influence  of  z, 
whereas  I*  is  in  the  zone  of  influence  of  w. 

(ii)  If  the  edge  e  lies  above  L,  Q  will  consist  only  of  I* 7m,  where  I*  is  defined  as 

above.  (Note  that  in  this  case  I*  is  not  entirely  contained  in  U;  however  by  allowing  I*  to 
intersect  the  exterior  of  U  in  such  cases,  we  simplify  the  handling  of  the  list  Q  without  losing 
accuracy  of  the  shortest  path  tracing  which  is  what  we  require.  A  similar  observation  applies 
to  cases  in  which  no  vertex  is  encountered,  but  in  which  the  upper  and/or  lower  boundaries 
of  U  progressively  invade  and  cut  off  more  and  more  intervals  of  the  list  Q.) 

Symmetric  rules  govern  the  case  where  z  lies  on  the  upper  boundary  c2  of  U. 

These  observations  allow  us  to  maintain  the  list  Q(t)  in  the  following  manner:  Intervals 
of  Q  can  be  associated  with  the  points  in  /V  in  whose  zone  of  influence  they  lie  and  their 
endpoints  represented  by  the  coefficients  of  the  rays  on  which  these  endpoints  lie.  The 
resulting  combinatorial  description  changes  only  when  /(f)  passes  through  one  of  the  corners 
of  fV,  which  we  can  suppose  to  have  been  presorted  according  to  their  *- coordinates  in  left  to 
right  order.  We  process  each  z  €  N  in  turn.  If  z  lies  on  cu  we  search  Q  in  linear  order  from 


-9- 

its  lowest  interval  upwards,  until  the  interval  /  containing  z  is  found.  Let  w  be  the  point  in  N 
in  whose  zone  of  influence  /  lies.  We  assign  v(z)  :=  w,  truncate  the  lower  portion  of  /  at  z, 
thereby  obtaining  a  new  /*,  and,  if  the  line  passing  through  w  and  z  has  a  slope  greater  than 
that  of  the  edge  of  cx  proceeding  forward  from  z,  add  a  new  interval  70  below  /*.  All 
intervals  of  Q  encountered  during  the  ascending  search  for  the  interval  containing  z  are 
deleted  from  Q.  Similar  and  symmetric  update  rules  are  applied  if  z  lies  on  c2. 

In  this  manner  we  obtain  the  values  v(z)  for  all  z  €  N.  In  processing  the  final  point 

i 

B  €  N  no  updating  of  Q  is  required;  we  only  need  to  scan  Q  to  find  the  interval  containing  B, 
and  hence  also  v(B). 

Once  the  map  v  becomes  available,  the  corners  of  the  required  shortest  path  p  are  (in 
reverse  order)  the  points  xx,  .  .  .  ,xd,  where  x1  =  B,  xd  —  A,  and  for  each  »  <  d  we  have 
x(+i  =  v(xt)-  This  path  is  easily  obtained  by  tracing  the  map  v  from  B  backwards  until  we 
reach  A . 

In  regard  to  the  running  time  of  this  algorithm,  note  that  if  an  interval  /'  contained  in 
the  zone  of  influence  of  some  z  6  N  is  deleted  from  Q(t),  then  the  zone  of  influence  of  z  will 
never  reappear  in  Q  for  any  t'  >  t,  and  thus  the  number  of  deletions  from  Q  is  bounded  by 
the  total  number  n  of  corners  in  N.  It  is  equally  plain  that  the  total  number  of  truncations  of 
intervals  in  Q ,  and  the  total  number  of  insertions  of  new  intervals  into  Q  are  both  bounded 
by  n.  Since  each  deletion,  truncation,  and  insertion  can  be  performed  in  constant  time,  it 
follows  that  the  shortest  path  algorithm  that  has  just  been  described  runs  in  time  linear  in  n. 

Remark:  An  interesting  application  of  this  smoothing  procedure  is  as  follows.  Suppose  that 
v  =  c*(x)  is  a  noise-perturbed  approximation  to  a  smooth  function,  and  that,  given  a 
sequence  of  points  lying  on  the  graph  of  c,  we  wish  to  find  the  regions  of  convexity  and  of 
concavity  of  c.  Apply  the  above  smoothing  algorithm  to  the  strip  U  bounded  between 
c:(x)  =  c*(x)  -  t  and  c2(x)  =  c*(x)  +  e,  for  sufficiently  small  t.  It  can  be  shown  that  if  the 
given  sequence  of  points  on  c  is  sufficiently  dense  then  good  approximations  to  the  regions  of 
convexity  can  be  characterized  as  regions  in  which  the  shortest  path  c'  within  U  between  the 
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Midpoints  of  c  passes  through  corners  in  the  lower  boundary  cx  of  U,  whereas  regions  of 
concavity  will  be  very  nearly  those  in  which  c'  passes  through  corners  of  c2.  This  method  is 
preferable  to  direct  comparison  of  the  chord  connecting  two  points  ylf  y2  (lying  sufficiently 
close  together  on  c)  to  their  midpoint  x  one,  because  c'  eliminates  much  of  the  noise  present 
in  the  raw  data. 

Smoothing  an  arbitrary  simple  polygonal  plane  curve. 

The  basic  algorithm  just  described  applies  only  to  cases  in  which  the  curve  c  being 
approximated  is  monotonic  in  at  least  one  coordinate  direction.  When  this  does  not  hold  for 
the  whole  of  c,  a  simple  way  of  proceeding  is  to  decompose  c  into  disjoint  sections,  each  of 
which  is  monotonic  in  some  direction,  and  then  to  apply  the  above-described  shortest  path 
algorithm  to  each  of  these  sections  separately.  In  the  experiments  reported  below,  we 
employed  this  simple  technique,  and  more  precisely  have  used  a  decomposition  scheme 
involving  four  directions  of  potential  monotonicity  -  the  x  axis,  the  y  axis,  and  the  two  lines 
at  45  degrees  to  these  coordinate  axes.  Given  a  curve  c,  we  find  the  longest  initial  portion  of 
c  which  is  monotonic  in  at  least  one  of  these  four  directions.  After  breaking  off  this  initial 
section,  we  find  the  longest  initial  portion  of  the  remainder  of  c  which  is  monotonic  in  one  of 
these  directions,  and  continue  in  this  manner  until  all  of  c  is  processed.  This  simple  heuristic 
proves  to  be  adequate  even  in  the  presence  of  relatively  large  amount  of  noise. 

Remark:  As  noted  above  there  exist  algorithms  due  to  Lee  and  Preparata  [LP84]  and  to 
Chazelle  [Ch82],  for  finding  the  shortest  path  between  two  specified  points  inside  any  given 
simple  plane  polygon.  These  algorithms  run  in  time  0(n  log  n),  where  n  is  the  number  of 
vertices  of  the  polygon.  A  difficulty  in  applying  these  algorithms  to  the  case  of  an  arbitrary 
simple  polygonal  curve  is  that  it  is  not  clear  what  is  the  best  way  to  define  the  polygon  within 
which  the  shortest  path  will  be  sought.  For  this  reason,  simple  extension  of  the  linear-time 
algorithm  given  above  for  monotonic  curves  seems  a  suitable  approach. 
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Approxlmatlon  by  nearly-shortest  paths  in  three  dimensions 

Next  we  generalize  the  preceding  considerations  to  curves  in  E3.  Since  in  this  case  no 
adequately  efficient  technique  for  finding  shortest  paths  lying  in  polyhedral  regions  is  known 
(cf.  [SS84b],  [Pa84]),  we  will  use  a  modified  shortening  technique  for  smoothing  which 
attains  the  efficiency  required.  As  before,  let  c  be  a  smooth  curve,  but  now  suppose  that  c  is 
three-dimensional,  and  let  c*  be  obtained  from  c  by  addition  of  a  perturbing  (noise)  function 
of  maximum  magnitude  at  most  e.  For  simplicity,  we  will  suppose  that  c  (resp.  c*)  has  a  1-1 
parametrization  of  the  form  (x,y(x),z(x))  (resp.  (x ,y* (x) ,z* (x))) .  Define  the  neighborhoods 
t/,(c)  and  t/t(c*)  in  terms  of  'Manhattan'  distance,  e.g.  Ut(c*)  is  the  set  of  all  points  (x,y0,z0) 
such  that  x  belongs  to  the  (common)  domain  0  ^  x  ^  X  of  c  and  c*,  while 

y*(x)  -cSji0S  y*(x)  +  e  , 

z*(x)  -  t  s  r0  s  z*(x)  +  e  . 

The  projection  of  Ut(c*)  into  the  (x.y)-plane  is  the  set  Vt(cm)  of  all  (x,y0)   satisfying 

0  s  x  =s  X  and  y'(x)  -«sy0s  y*(x)  +  e.    Define  the  function  y:(x)  for  0  ^  x  •£■  X  by 

requiring  that  the  curve  (x,V;(x))   shall  be  the  shortest  curve  in  Vt{c*)  connecting  the 

projections  onto  the  (x,v) -plane  A,  B  of  the  endpoints  a,  b  of  c*.  Since  V€(c')  is  a  2-D  'band' 

of  the  form  considered  above,  the  planar  shortest  path  algorithm  that  we  have  described 

calculates  y^x)  efficiently. 

Next  let  W^c'.yJ  be  the  set  of  all  points  (x,y1(x),z0)  such  that  0  ^  x  ^  X  and 
2*(x)  -  t  s  z0  <  z*(x)  +  €,  and  define  the  function  z2(x)  for  0  =s  x  s  X  by  requiring  that 
the  curve  (x,y1(x),r1(x))  shall  be  the  shortest  curve  in  Wt(c*,yJ  connecting  the  endpoints  of 
c*.  Since  Wt(c*,y{)  is  not  too  different  from  the  plane  (it  is  in  fact  a  2-dimensional  piecewise 
planar  surface  folded  at  those  x  which  are  the  corners  of  y^x))  our  planar  shortest  path 
algorithm  also  suffices  to  calculate  zx(x)  efficiently  (simply  apply  this  algorithm  to  the  planar 
band  obtained  by  unfolding  Vt(cm,yJ).  We  shall  show  that,  provided  €  is  small  enough,  the 
curve  c^x)  =  (x  ,y  :(x)  ,z  :(x))  is  a  good  approximation  to  the  original  curve 
c(x)  =  (x,v(x),z(x)),  in  the  sense  that  after  both  these  curves  are  parametrized  by  arc-length 
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they  must  approximate  each  other  uniformly. 

For  this  we  can  employ  much  the  same  argument  used  to  prove  Lemma  2.2,  provided 
that  we  first  show  that  the  curve  cl  is  not  much  longer  than  the  model  curve  c.  Note  that 
since  cx  lies  within  U2t(c),  the  length  of  cx  must  be  at  least  L  -  8,  since  we  mean  to  assume 
that  c  has  length  L  and  is  (2t,5)-stable. 

To  prove  that  cx  is  not  much  longer  than  c,  we  argue  as  follows.  In  a  neighborhood  O 
of  the  projected  planar  curve  C  =  (x,y(x))  introduce  local  curvilinear  coordinates  (s,u)  which 
map  C  onto  the  x-axis,  s  being  arc-length  along  the  curve  C,  and  which  map  each  sufficiently 
short  segment  orthogonal  to  the  curve  C  into  a  segment  of  the  form  (s0,u), 
-a(sQ)  s  u  ^  b(s0).  (Here  it  is  plain  that  a  and  b  must  be  positive  functions.)  In  these 
curvilinear  coordinates  the  Euclidean  (.r,y) -length  of  any  plane  curve  p(t)  is  expressed  by  an 
integral 

J<p'(0gCp(0)i»'(0t)1'2*.  (i) 

where  as  usual  p'  (r)  designates  the  derivative  of  the  (vector- valued)  function  p(t),  and  where 
G{p)  is  a  positive  definite  2x2  matrix.  Moreover,  since  arc-length  along  the  j-axis  is 
measured  directly  by  the  coordinate  s,  and  since  the  lines  (s0,u)  are  all  orthogonal  to  the  .i- 
axis  in  the  metric  defined  by  G,  we  must  have  G((s,0))  -  I  identically  in  s,  where  / 
designates  the  2  x  2  identity  matrix.  It  follows  that  for  each  positive  constant  -q  we  have 

v  G((s,u))  vT  a  (l--n)2v  vr 
for  all  vectors  v  and  all  sufficiently  small  u. 

If  e  is  sufficiently  small,  the  minimum  length  curve  (jt.y^jt))  in  Vt(c*)  introduced  above 
can  be  written  in  (s,u)  coordinates  as  (s,y(.s))  where,  unfortunately,  the  function  Y  may  be 
multi-valued.  (However,  by  perturbing  y^.r)  very  slightly,  we  can  assume  that  only  a  finite 
number  of  smooth  'branches'  of  Y{s)  lie  over  each  interval  of  the  j-axis.)  This  will  cause 
slight  technical  complications  but  will  be  seen  not  to  disturb  the  overall  course  of  the 
argument  we  are  about  to  give.  Let  Lc  be  the  length  of  the  curve  C.  Since  C  is  included  in 
the  projected  neighborhood  Vt(c*)  of  c*,  the  length  of  the  curve  (x.y^jc))  is  no  more  than 
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Lc.  Expressing  this  fact  in  the  curvilinear  coordinates  (s,u),  we  find  that 

/  (  (l,rW)  G((s,Y(s)))  (l,r(*))r)W  ^Lc,  (2a) 

o 

Hence,  assuming  that  e  is  sufficiently  small  relative  to  the  positive  constant  t\  chosen  above, 
we  have 

(l-^/fl  +  rWfAsLf,  (2b) 

0 

i.e. 

/  [  (1  +  r(*)2)w  -  1  ]  *  s  —3—  Lc  .  (2b) 

0  X   —     M 

Note  that  the  inequalities  (2a)  and  (2b)  hold  exactly  as  written  if  the  function  Y(s)  is  single 
valued,  whereas  if  Y(s)  is  multi  valued  and  represented  by  multiple  branches  over  certain 
intervals  of  the  j-axis,  the  (positive)  integral  in  (2a)  and  the  integral  in  (2b)  should  be 
understood  to  represent  the  sum  of  corresponding  single  valued  integrals,  the  sum  being 
extended  over  all  the  branches  of  Y(s). 

It  is  easy  to  see  that  the  following  inequalities  hold: 


(1  +  a2)1'2  -Uja3,       0  s  a  s  1  , 
(1  +  a2)1'2  -Uja,       lSo. 


It  therefore  follows  from  (2b)  that 

/  |r(j)P*s3-r2L_Lc 

irtJiisi  *  -  "n 


(3a) 


and 


/  \Y'(s)\ds*3-^—Lc.  (3b) 

As  before,  these  inequalities  are  valid  exactly  as  written  if  the  function  Y(s)  is  single  valued, 
whereas  if  Y(s)  is  multi  valued  and  represented  by  multiple  branches  over  certain  intervals  of 
the  5-axis,  each  of  the  integrals  in  (3a)  and  (3b)  should  be  understood  to  represent  the  sum 
of  corresponding  single  valued  integrals,  the  sum  being  extended  over  appropriate  portions  of 
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aU  the  branches  of  Y(s). 

Our  next  step  is  to  show  that  the  curve  c3(x)  =  (x,yi(x),z(x)),  0  s  j  s  x,  is  not  much 
longer  than  the  curve  c(x)  =  (x,y(x),z(x)),  i.e.  has  length  not  much  greater  than  L.  To  see 
this,  parametrize  c3  in  terms  of  s  as  (s,Y(s),Z(s)),  where  just  as  above  the  function  Y  is 
multi-valued  and  can  therefore  cause  slight  technical  complications.  (However,  Z(s)  is 
single-valued  since  it  merely  represents  z(x)  which  is  single-valued  along  the  curve  c,  i.e.  the 
j-axis.) 

The  positive-definite  matrix  G  appearing  in  (1)  satisfies  G((s,u))  s  (1  +  t|)  /  for  all 
sufficiently  small  u.  Thus,  assuming  as  above  that  c  is  sufficientiy  small,  we  have 

/  (  (l,rW)  G«s,Y(s)))  (l,Y'(s))T  +  Z'(s)2  Y2  ds  (4) 

o 

s(i  +  1)/(i  +  r(,f+z'(i)2)^) 

0 

where  as  above  we  may  be  required  to  sum  in  appropriate  fashion  over  multiple  branches  of 
Y(s).  The  right  hand  side  of  the  inequality  (4)  can  be  written  as 


[<x                                             ^c 

(1  +  T|) 

J(1  +  Z'(,):)L'2*  +  J 

0                                                  0 

the  second  in 

tegral  B  in  (5)  satisfies 

(i  +    r(*f  y2  _  i 

1  +  Z'(s)2 


(1  +  Z'(s)2)v2ds 


(5) 


**j       J  \Y'(s)\2ds+       f  \Y'(s)\ds^5-^—L,  (6) 

z   |J"(i)|s  1  |l"(f)|>  1  i  -  ^1 

by  (3a)  and  (3b),  since  (1  +  a2)1'2  -  1  <  y  a2  and  (1  +  a2)1-'2  -  1  <  a  for  all  a  ;>  0.    As 

before,  let  L  be  the  length  of  the  model  curve  c  =  (x,  y(x),  z(x)),  so  that  plainly  Lc  s  L 
Since  the  first  integral  in  (5)  (resp.  the  integral  in  the  left-hand  side  of  (4))  is  simply  Lc 
(resp.  the  length  L3  of  c3,  we  find  that 

L3  ^  (1  +  „)  L  +  5  iO+Jll  L  ,  (7) 

1  -  tj 

showing,  as  required,  that  when  t  is  sufficiently  small,  L3  is  not  much  larger  than  L. 

Since  the  curve  c3(x)  =  (x,y1(x),z(x))  belongs  to  the  set  Wt(c*,yJ  introduced  above, 
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and  since  c:(x)  =  (x.v^jO.z^jc))  is  the  shortest  curve  in  Wt(c*,yJ  connecting  endpoints 
which  differ  only  slightly  from  those  of  c3,  it  follows  that  ci  also  has  length  differing  only 
slightly  from  that  of  c.  Thus,  as  noted  above,  c:  will  be  uniformly  good  approximation  to  c  if 
€  is  sufficiently  small  and  both  c  and  c:  are  parametrized  by  their  arc  length. 

The  preceding  discussion  assumes  that  the  whole  of  the  given  curve  c  can  be 
parametrized  by  the  coordinate  x.  To  treat  more  general  simple  3-D  curves  we  can 
decompose  therq  into  raonotonic  subsections  using  much  the  same  heuristic  described  for  2-D 
curves,  and  then  apply  the  shortening  procedure  just  described  to  each  one  of  these  sections. 


Fig.  2.1  and  Fig.  2.2  show  a  2-D  curve  artificially  perturbed  by  addition  of  random 
noise,  and  the  result  of  smoothing  this  curve  by  the  'shortening*  method  just  described;  figs. 
4.1  and  4.2  of  section  4  do  the  same  for  a  3-D  curve. 


Fig.  2.1.  An  Artificially  Perturbed  2-D  Curve. 
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Fig.  2.2.  Fig.  2.1  Curve  After  Smoothing. 
3.   Fast  boundary  matching  by  least  squares  fitting. 

Next  we  describe  a  fast  technique  for  matching  a  portion  of  the  boundary  C"  of  an 
observed  image  O'  against  a  corresponding  portion  of  the  boundary  C  of  a  model  object  0. 
We  assume  that  C"  has  been  smoothed  using  one  of  the  smoothing  methods  described  in  the 
preceding  section,  and  that  both  C  and  C"  are  parametrized  by  arc  length.  Recall  that  the 
matching  we  seek  calls  for  determination  of  the  offset  s0  and  the  Euclidean  transformation  E 
for  which  the  curves  EC(s  +  s0)  and  C  (s)  are  closest  to  one  another  in  the  L2  norm.  To  be 
more  specific,  represent  each  of  the  curves  C,  C  by  a  sequence  of  evenly  spaced  points  on  it, 
and  let  these  sequences  be  (<iy)"-i  and  (vy)"_i,  corresponding  to  an  observed  (smoothed) 
object  0'  and  a  model  object  O  respectively.  For  simplicity  in  describing  the  necessary 
calculation,  assume  first  that  the  whole  boundary  of  O'  is  visible  and  that  no  offset 
calculation  is  required.  Matching  then  amounts  to  finding  the  Euclidean  motion  £  of  the 
plane  which  minimizes  the  I2  distance  between  the  sequences  (Eu,)jm[  and  (v,)?.^  i.e.  we 
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need  to  compute 

n 

A  =  min  2  [Euy  -  vy|2 
To  simplify  the  calculation,  we  first  translate  O'  so  that  the  centroid  of  its  boundary  lies  at 
the  origin,  i.e.  so  thai 

Next  we  write  E  as  En  =  R6u  +  a,  R9  denoting  a  counterclockwise  rotation  by  9.  Then 


A  =  min    2  [/?9uy  +  a  -  vy|2  = 


«.»      y-i 


min 


2  |vy|2    +  n|a|2  -  22  a-v,  +  2  l«yl2  +  2£  «*e°y  "  2£  *.uyvy 

L/-i  y-i  y-i  y-i  y-i 


But 


2  a-/?9uy  =  a/?9(2uy)  =  0  . 
Hence  a  and  8  appear  independently  in  A  and  we  can  minimize  their  contributions  separately. 


To  minimize  over  a  we  simply  put 


1     " 


As  to  0,  we  need  to  compute 


8  =  max   2  ^eu/vy 

9     y-i 

Regarding  the  vectors  uy,  vy  as  complex  numbers  Uj,  Vj,  we  can  rewrite  this  as 


8  =  max  Re 
e 


2  *'V, 


=  IS  ujvj\  > 


y-i 


where  plainly  the  maximizing  9  is  just  the  negative  of  the  polar  angle  of  2  ujvj-    Altogether 
this  gives 


A  =  2  N2  "  ht  v,P  +  ±  |nyP  -  2|£  «yvy| 


C) 


y-i 


y-i 


y-i 


y-i 


A  =  2  |v,|2  -  1|£  vy|2  +  ±  |u/  -  2  |  £  uy-vy  |2  +  |  ±  uyxvy  P 
y-i  n  y-i  y-i  I  y-i  y-i 

where  uxv  denotes  the  (2-dimensional)  cross  product  of  the  vectors  u  and  v. 
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If  O'   is  partially  occluded,  we  have  to  match  the  sequence  (uy)"_i  to  each  of  the 
contiguous  subsequences  (▼/+(f)"-i  of  the  (circular)  sequence  (▼/)£.  i,  for  d  =  0,  .  .  .  ,m-\. 
(Here  we  can  appropriately  assume  that  m  a  n,  for  otherwise  the  (partial)  periphery  of  O'  is 
too  long  to  match  O.) 
For  each  such  d  (*)  becomes 

aw  =  2  N2  -  7I  2  v,l2  +  2  N2  -  2I2  «,v„| 

It  is  easily  seen  that  the  minimum  of  the  values  A(«f),  rf  =  0,  .  .  .  ,m-l,  can  be  found  in 
time    O(mlogm),    by    using    the    fast   Fourier    transform    algorithm    for   computing    the 

n 

convolutions  2  ujVj+d- 
/-i 

The  matching  technique  just  described  generalizes  easily  to  curves  in  three  dimensions, 
or,  more  generally,  to  any  situation  in  which  a  model  curve  or  surface  u(io)  depending  on 
one  or  more  parameters  id  must  be  rotated  and  translated  to  match  a  model  curve  or  surface 
v(cj)  as  well  as  possible.  We  must  assume,  however,  that  the  matching  operation  involves  no 
change  in  parametrization  for  either  of  the  functions  u(a>)  or  v(a)). 

Suppose  then  that  we  are  given  two  descriptor  functions  u(u),  v(cj),  w  £  S, 
corresponding  respectively  to  an  observed  curve  or  surface  C"  and  a  model  curve  or  surface 
C.  We  need  to  find  the  Euclidean  motion  E  (e.g.  of  3-space)  which  minimizes 

A  =  min   /  |£ u(a>)  -  v(u>)|:  dtx> 
E      s 

(The  domain  S  of  integration  can  have  one  or  more  dimensions.)  As  in  the  2-D  case,  we 
translate  C"  so  that  its  centroid  lies  at  the  origin,  giving 

J  u(u>)  dm  =  0 

5 

Write  E  as  £u  =  /?u  +  a,  where  R  is  a  rotation.  Then 

A  =  min    J  \Ru(w)  +  a  -  v(w)p  du  = 

/  |v(<i))p  doi    +  \S\\a\2  -  2  /  a-v(to)  dta  +  J  |u(co)|2  dot  +  2/  a-/?u(o>)  du>  -  2  f  /?u(u>)v(u>)  dm 
s  s  s  s  s 
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But 

/  a  Ru((i})  du)  =  *R(J  u(a>)  du  )  =  0  . 

5  5 

Hence   ■   and  R    appear   independently   in  A   and  we  can  minimize  their  contributions 
separately. 

To  minimize  over  a  we  set 


As  to  R ,  we  need  to  compute 


■  -  "TjT  /  *(o)  du 


8  =  max    f  Ru(oi)\(oi)  du 

R       s 

To  find  8  we  first  calculate  the  matrix  A  given  by 

Atj  -  /  U/(">)Vy(o>)  du  . 

5 

(where  ij  =  1,2,3  if  we  are  dealing  with  a  curve  or  surface  in  3-space).  In  terms  of  the 
matrix  A  we  can  express  S  as 

8  =  max    tr(RA) 

R 

To  maximize  tr(RA)   first  assume  A    to  be  non-singular;  then  we  can  decompose  A   as 

A  =  QH,  where  Q  =  A(A*A)  2  is  a  pure  rotation  and  H  =  (A* A)2  is  positive  definite 
symmetric.   This  gives 

i 
8  =  max   tr(RQH)  -  max   tr(RH)  =  tr{H)  =  tr((A*A)2) 

X  A 

Indeed,  since  the  trace  is  invariant  under  rotation,  we  may  assume  that  H  is  diagonal.  But 
since  for  a  diagonal  positive  definite  matrix  (X,)  and  a  rotation  matrix  (r,j)  the  trace  of  their 
product  2  ^tru  can  De  no  larger  than  5)  *;  and  can  assume  this  value  only  when  r»  =  8y. 
An  identical  result  for  a  singular  A  follows  easily  by  continuity  arguments.  Altogether  we 
have 

,  I 

A  =  /  |v(«)P  du  -  -i-IJ  v(co)  rfwp  +  /  |u(o»)|2  du  -  2  tr((A'A)2)  (") 

S  PI     5  S 
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showing  that  the  optimal  rotated  match  between  C  and  C  can  be  found  in  time  proportional 
to  that  needed  to  integrate  the  various  functions  appearing  in  (••),  i.e.  proportional  to  the 
number  of  data  points  used  to  discretize  the  curves  or  surfaces  o  and  v.  In  particular,  as  in 
the  2-D  case  considered  previously,  these  formulae  can  be  used  to  match  observed  3-D  curves 
parametrized  by  arc  length  to  similarly  parametrized  model  curves.  Matching  can  be  achieved 
in  time  0(n  log  /»)  by  using  the  fast  Fourier  transform,  even  if  the  observed  curve  O'  is 
partially  obscured.  Formula  (**)  gives  us  the  smallest  least-squares  distance  between  the 
objects  C  and  C",  but  does  not  produce  the  Euclidean  transformation  E  which  attains  this 
best  match.  In  case  A  is  nonsingular,  the  required  E  is  given  by  EU  =  R  u  +  a  where  a  is  as 
given  above,  and  R  -  (A"  A)1'2  A~l.  Even  if  A  is  singular,  it  can  be  factored  as  A  =  QH 
where  Q  is  a  rotation  and  H  is  as  above,  but  in  this  case  any  rotation  of  the  form 
R  =  ^ofi"1'  where  R0  is  the  identity  on  the  range  of  H  but  can  be  arbitrary  on  the  space 
orthogonal  to  this  range,  will  minimize  tr(RA).  The  map  Q  can  be  constructed  in  this 
singular  case  by  finding  a  maximal  linearly  independent  set  of  columns  of  H,  and  then 
forming  Q  which  maps  each  such  column  into  the  corresponding  column  of  A  and  which  is 
the  identity  on  each  element  of  the  (common)  nullspace  of  A  and  H. 

4.   Application  to  the  analysis  of  compound  2-D  scenes 

A  'compound'  2-D  scene  is  one  consisting  of  thin  objects,  e.g.  objects  cut  out  of  metal 
or  cardboard,  which  may  overlap.  To  apply  a  boundary  matching  procedure  like  that 
described  in  the  preceding  section  to  such  scenes,  one  needs  a  way  of  determining  where  the 
boundary  of  one  of  several  overlapping  objects  ends  and  that  of  another  begins.  For  obvious 
reasons,  these  'breakpoints'  are  likely  to  be  points  at  which  the  tangent  to  the  overall 
boundary  B  of  the  overlapping  objects  changes  rapidly,  hence  points  at  which  the  second 
derivative  of  the  boundary  curve  reaches  a  peak.  To  be  more  specific,  let  b(s)  be  a 
parametrization  of  the  (smooth)  boundary  B  by  arc  length,  and  suppose  that  b(s)  has  the 
additional  property  that  as  we  traverse  B  in  the  direction  of  increasing  s,  the  interior  of  the 
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( union  of  the)  overlapping  objects  lies  to  the  right  of  B.  Let  n(s)  be  the  outward-pointing 
unit  normal  to  B  at  b(s).  Then  the  second  derivative  b"(s)  is  parallel  to  n(s),  and  we  look 
for  points  s  where  i2C0  ™  b"(s)  '  «(*)  reaches  a  positive  peak.  (Only  positive  peaks  are  of 
concern  since  negative  peaks  correspond  to  sharp  boundary  convexities  rather  than 
concavities,  which  are  not  likely  to  be  points  where  the  boundaries  of  two  distinct 
overlapping  objects  meet.) 

Specifically,  assume  that  the  boundary  B  has  been  smoothed  and  that  we  have 
enumerated  a  sequence  of  evenly  spaced  points  in  some  plausible  clockwise  order  around  5. 
Then  we  can  proceed  in  the  following  manner. 

(a)  Choose  some  fixed  (and  relatively  small)  k  in  advance,  take  it  successive  boundary  points 

starting  at  each  successive  place  j,  and  estimate  the  boundary  tangent  by  calculating  the  line 

of  best  (least  squares)  fit  to  these  successive  points. 

■ 

(b)  Calculate  the  second  derivative  by  differencing  successive  tangents   (see  below  for 

details),  and  look  for  maxima  of  the  second  derivative.  (As  noted,  by  ignoring  minima  we 
bypass  convex  corners  and  record  only  concave  corners.) 

(c)  If  many  maxima  lie  close  to  each  other,  thin  them  out. 

The  necessary  calculations  are  fast  and  simple.  Let  the  k  points  presently  being 
processed  be  p,  =  (xhy,),  i  =  1,  .  .  .  ,k.  The  horizontal  line  of  best  approximation  is  at  the 
level  v  which  minimizes 


j  2  (y/  -  y)2 > 

hence  v  =  —  S  v/>  anc^ tne  minimum  is 


*    i 


}Xv,2-(}2v;)2. 

This  analysis  applies  to  lines  of  any  orientation,  and  thus  the  line  of  best  approximation  must 
pass  through  the  center  of  gravity  of  the  set  of  points  p,,  and  the  unit  normal  v  to  this  line 
must  be  that  which  minimizes 
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H(v)  =  k  2  (ypt?  -  (v  2  P,)2  =  v  (P  -  Q  o  Q)  v  , 
where  P  is  the  dyadic  sum  2  P;  op„  and  G  =  2  ^/-  The  minimizing  vector  is  therefore  the 

eigenvector    belonging    to    the    smaller    eigenvalue    of    the    symmetric    2x2    matrix 

A  =  P  —  Q  oQ,  which  we  can  write  as 


(::)• 


whose  characteristic  equation  is  X2  -  (a  +  c)\  +ac  -  b2,  giving  the  smaller  eigenvalue 

a  +  c  -   [(a  -  c)2  +  4b2) 
X  =  2 

and  a  corresponding  eigenvector  proportional  to  both  (—b,  a  —  \)  and  (c-X,  —b).  (Care  has 

to  be  taken  to  use  the  larger  of  these  two  vectors,  so  as  to  avoid  situations  in  which  the 

smaller  of  the  two  vanishes.) 

The  calculations  required  can  therefore  be  pertormed  quite  efficiently  by  a  linear  scan 
through  the  sequence  of  boundary  points. 

To  find  concave  corners  using  these  best-fitting  tangents  we  can  take  the  cross  product 
of  each  such  vector  t  with  the  next  tangent,  which  gives  an  estimate  for  the  angle  between 
them.  Then  we  can  look  for  positive  peaks  in  this  angle  at  which  it  exceeds  some  threshold 

angle  — ,  the  denominator  k  being  appropriate  since  least  square  fitting  averages  the  tangent 

slope  among  k  successive  points,  making  it  reasonable  to  expect  that  roughly  *  steps  will 
transpire  before  one  sees  all  of  any  fast  change  in  tangent  orientation. 

Simulations  in  2  and  3  dimensions;  Experiments  with  real  2-D  images 

This  technique  just  explained  has  been  simulated  in  combination  with  the  smoothing 
method  described  above.  It  appea:s  to  work  well  and  to  withstand  considerable  noise  in  the 
boundary  curves  simulated.  Figures  4.1  and  4.2  show  the  result  of  one  such  simulation.  Fig  . 
4.1  figure  shows  (a  3-D)  a  test  curve  and  the  result  of  perturbing  a  translated  and  rotated 
subsection  of  it  by  artificial  noise.  Fig.  4.2  shows  the  matching  of  the  noisy  figure  to  the 
original  curve  obtained  after  smoothing  using  the  3-D  smoothing  procedure  described  above. 
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Fig.  4.1 .  A  Test  Curve  And  A  Perturbed,  Rotated  And  Translated  Subsection. 


The  companion  paper  [KSSS85]  reports  on  various  much  more  extensive  experiments  with 
partially  obscured  2-D  curves  and  composites  of  such  curves,  which  show  that  the  techniques 
described  above  can  be  used  to  analyze  such  configurations  robustly  and  reliably. 
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Fig.  4.2.   Fig.  4.1  Curve  Smoothed  And  Matched  To  Model  Curve. 


5.   Characteristic  Curves  on  the  Visible  Surfaces  of  3-D  Bodies 


To  apply  the  matching  method  described  in  the  preceding  pages  to  partially  obscured  3- 
D  bodies  0,  the  most  desirable  approach  is  to  define  and  match  suitable  characteristic  curves 
on  their  surfaces.  (If  such  curves  are  not  available,  we  will  have  to  match  observed  to  model 
2-D  surfaces,  which  may  be  a  much  more  expensive  operation,  especially  if  we  do  not  know 
which  portion  of  a  3-D  model  object  has  been  observed.)  Surface  curves  used  for  recognition 
must  be  curves  C  defined  by  the  geometry  or  surface  color  of  a  surface  O  and  fixed  in  O,  so 
that  when  O  undergoes  a  rotation  and/or  translation  the  curve  C  moves  along  with  it.  C's 
definition  should  depend  only  on  local  geometric  and  other  features  of  O,  so  that  parts  of  C 
can  be  constructed  even  when  only  a  portion  of  0's  surface  can  be  observed.  Such 
characteristic  curves  can  be  defined  in  any  one  of  the  following  ways: 

(i)  Use  of  albedo  or  color  variations.  Suppose  that  the  surface  of  O  is  'painted',  e.g.  that 
portions  of  this  surface  have  sensibly  different  reflectances  in  regions  of  the  electromagnetic 
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spcctrum.  A  simple  example  might  be  a  white  sphere  or  other  convex  body  with  red  curves 
painted  on  it.  Then  any  sharply  defined  reflectance  transitions  can  be  used  as  characteristic 
curves.  Moreover,  even  if  no  sharp  transitions  exist  because  O's  reflectance  varies  smoothly 
over  its  surface,  it  may  still  be  possible  to  use  the  'level'  curves  along  which  this  reflectance 
has  a  given  value  as  characteristic  curves. 

Note  that  we  assume  here  that  a  3-D  or  'depth'  image  of  the  visible  surfaces  of  O  is 
available,  so  that  the  true  3-D  geometry  of  each  characteristic  curve  C  observed  on  O's 
surface  is  known.  Such  information  can  easily  be  obtained  by  combining  any  suitable  depth 
sensor  with  a  standard  video  camera  sensitive  to  appropriate  parts  of  the  optical  spectrum. 
(When  color  variations  on  O's  surface  are  being  used,  it  may  be  advantageous  to  eliminate 
Lambert-law  dependencies  of  reflectance  on  the  angles  of  illumination  and  observation  by 
forming  ratios  of  two  or  more  separate  images  acquired  using  light  of  various  wavelengths.) 

(ii)  'Curves  of  rolling'.  Next  consider  the  somewhat  more  difficult  case  in  which  O's 
surface  is  'pure  white',  i.e.  in  which  only  geometric  information,  but  no  useful  reflectance 
information,  is  available  concerning  these  surfaces.  Even  in  this  case  it  is  not  hard  to  define 
characteristic  curves.  The  easiest  case  to  consider  is  that  in  which  0  is  not  convex,  in  which 
case  'curves  of  rolling'  on  O's  surface  will  exist  and  can  be  used  as  characteristics.  These 
curves  are  defined  as  follows.  Each  tangent  plane  to  a  smooth  convex  object  O  touches  it  in 
just  one  point,  but  if  0  is  not  convex  it  will  have  tangent  planes  P  which  touch  it  in  at  least 
two  separated  points.  (We  ignore  the  even  more  favorable  case  of  tangent  planes  which 
touch  O's  surface  in  three  separated  points,  since  in  general  there  can  exist  at  most  a  finite 
number  of  such  planes,  and  hence  as  soon  as  one  such  plane  has  been  found  the  problem  of 
identifying  0  and  determining  its  orientation  becomes  purely  discrete  (This  remark  applies  to 
polyhedra,  and  makes  them  easy  to  identify  by  a  relatively  straightforward  'probing'  method 
see  [GLP84],  [GLP85].) 

We  shall  say  that  a  point  p  on  the  surface  of  a  smooth  object  O  is  a  point  of  rolling  if  the 
tangent  plane  to  0  at  p  is  also  tangent  to  O  at  another  point  q,  and  call  the  locus  of  all  such 
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points  p  the  curves  of  rolling  on  the  surface  of  0.  (This  name  is  suggested  by  the  fact  that, 
given  a  physical  model  of  O,  curves  of  rolling  on  it  can  be  found  by  putting  O  on  a  flat 
tabletop  T  and  finding  positions  in  which  O  rolls  on  T  with  two  points  in  contact  with  T.  The 
points  at  which  O  contacts  T  as  it  rolls  along  the  one-parameter  path  available  to  it  are  then 
curves  of  rolling.) 

To  see  that  the  locus  of  points  of  rolling  on  a  smooth  3-D  body  0  can  generally  be 

expected  to  lie  along  curves,  we  can  argue  as  follows.    Map  each  point  p  of  O's  surface  S 

i 

into  the  plane  P(p)  tangent  to  O  at  p.  This  maps  the  two-dimensional  surface  S' ,  into  the 
space  II  of  all  planes  in  E3.  A  plane  is  delined  by  its  unit  normal  vector  and  distance  from  the 
origin,  three  parameters  in  all,  so  that  11  is  a  3-dimensional  manifold  into  which  P  maps  S. 
The  points  of  rolling  correspond  to  points  of  self-intersection  of  the  image  surface  P(S),  and 
thus  since  self-intersections  of  a  2-dimensionaJ  surface  embedded  (in  general  position)  in  a 
3-dimensional  space  like  II  must  lie  along  1-oimensional  curves,  it  follows  that  the  points  of 
rolling  on  O's  surface  will  generally  lie  along  smooth  curves. 

Note  that  points  of  rolling  can  be  found  as  soon  as  enough  of  O's  surface  to  show  the 
presence  of  two  convexities  separated  by  a  concavity  has  been  observed. 

(iii)  Other  geometrically  defined  invariants.  Even  if  0  is  convex  (or  if  its  concavities 
are  not  observed),  other  geometric  invariants  of  its  surface  can  be  used  to  define 
characteristic  curves.  Generally  speaking,  any  function  f(p)  defined  for  points  p  on  the 
surface  of  0  and  rotating  with  O  will  define  characteristic  'level  curves'  f(p)  =  it.  Another 
possibility  is  to  use  a  function  of  this  kind  which  is  defined  on  some  auxiliary  surface  O' 
associated  with  O  that  rotates  with  0.  For  example,  we  can  let  O'  be  the  sphere  of  unit 
vectors  v,  and  for  each  d  >  0  and  v  £  O'  define  fd{v)  as  follows:  take  a  plane  perpendicular 
to  v,  initially  very  far  out  in  the  direction  of  v,  and  then  translate  this  plane  P  back  along  v 
(maintaining  its  orientation)  until  it  first  touches  O  and  then  penetrates  into  O  to  depth  d. 
L^  fAv)  De  some  rotational  invariant  of  the  intersection  of  P  with  O,  e.g.  the  area,  central 
moment,  or  perimeter  of  this  intersection.    Another  possibility  (but  one  more  apt  to  be 
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sensitivc  to  noise)  is  to  map  each  point  p  of  O  to  its  normal  vector  N{p),  and  to  define  the 
characteristic  function  fd  by  some  approximation  to  the  Jacobian  of  the  map  N  (e.g.  the  area 
on  the  unit  sphere  swept  out  by  N{q)  as  q  traverses  the  part  of  O's  surface  which  lies  within 
distance  d  of  p). 

6.  A  surface-based  matching  procedure  for  recognition  of  a  general  class  of  'machined* 
objects. 

Having  discussed  curve-based  matching  procedures  in  the  preceding  sections,  we  now  go 
on  to  consider  what  can  be  done  by  way  of  matching  surfaces  on  which  stably  defined  curves 
cannot  easily  be  found.  Here  we  will  only  be  able  to  deal  efficiently  with  relatively  simple 
surfaces.  Fortunately,  however,  the  surfaces  of  machined  objects  are  often  quadric,  or,  if 
this  is  not  the  case,  can  be  represented  well  by  quadratic  or  cubic  splines.  The  following 
paragraphs  will  describe  a  matching  technique  for  such  surfaces,  which  can  be  made 
particularly  effective  for  quadric  surfaces.  The  technique  in  question  is  applicable  to  any 
surface  having  an  implicit  polynomial  representation  P(u,v,w)  =  0. 

Suppose  that  a  finite  collection  of  model  surfaces  with  such  descriptions  is  given,  that  a 
collection  of  points  x,  =  («,,  v;,  z,  )  on  one  of  these  surfaces  is  observed,  but  that  the 
Euclidean  orientation  E  of  the  observed  surface  is  unknown.  We  can  define  the  'degree  of 
match'  between  the  set  of  observations  and  the  equation  P(x)  =  0  for  the  model  as  the 
minimized  difference 

mjn  Y,wt[p(Ex,)Y,  (1) 

where  E  ranges  over  the  group  of  Euclidean  transformations,  and  where  the  wt  are  weights 
chosen  to  reflect  some  heuristic  notion  of  the  'evidential  weight*  of  each  observation.  (For 
example,  wt  can  be  the  distance,  or  the  squared  distance,  of  the  observed  point  xi  from  the 
average  of  all  observations;  this  will  increase  the  'evidential  weight'  assigned  to  outlying 
points,  which  may  be  appropriate). 

A  difficulty  is  that  (1)  involves  a  nonlinear  minimization  leading  to  equations  which  will 
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gencrally  be  of  unpleasantly  high  degree.  To  sidestep  this  difficulty,  we  can  minimize,  not 
over  the  map  £,  but  over  the  coefficients  of  all  polynomials  of  the  same  degree  as  the  model 
surface  P(x)=0  we  are  trying  to  match.  That  is,  we  can  find  the  polynomial  Q,  having  the 
same  degree  as  P,  at  which  the  minimum 

mm  2"i(e(*i))2.  (!') 

is  attained.   This  is  a  much  more  tractable  quadratic  minimization  problem.  Moreover,  if  the 

observed  points  did  lie  exactly  on  the  surface  P(Ex)=0,  then  the  minimizing  Q  in  (1*)  would 

clearly  be  Q(x)  =  P(Ex),  so  that  once  having  the  coefficients  of  Q  we  could  expect  to  calculate 

E  in  straightforward  fashion.    (Suppose,  for  example,  that  the  polynomial  P  is  of  even  order. 

The  rotational  part  R  of  E  acts  independently  of  £'s  translational  part  on  the  highest  order 

terras  of  P.    If  P  is  quadratic,  R  can  be  determined  from  Q  by  finding  the  principal  axes  of 

(2's  purely  quadratic  part  and  matching  them  with  the  corresponding  axes  for  P.  If  P  is  of 

higher  degree,  we  can  apply  the  Laplacian  operator  A,  which  is  rotationally  invariant,  to  the 

highest  order  terras  of  both  P  and  Q  just  often  enough  to  produce  two  quadratic  polynomials, 

and  then  match  principal  axes  as  before.    Next  suppose  that  the  principal  terms  of  P  are 

cubic.    In  this  case,  we  can  associate  an  orthogonal  set  of  'principal  axes'  with  P  in  the 

following  way.    Take  AP,  which  is  linear,  and  take  the  vector  v;  orthogonal  to  the  plane 

ikP-0.  This  is  the  first  principal  axis  of  P.  Then  take  the  (cubic)  restriction  of  P  to  the  plane 

AP  =  0,  and  repeat  this  procedure  recursively  to  define  n-1  additional  axes  v;,...,vn  which 

are  orthogonal  to  v:.  These  axes  have  an  invariant  relationship  to  P,  in  the  sense  that  they 

simply  rotate  when  P  is  transformed  into  P(Rx).    If  Q  and  P  are  both  cubic,  they  can  be 

matched  simply  by  matching  their  characteristic  axes;  if  they  are  of  higher  degree,  we  can 

apply  the  Laplacian  operator  an  appropriate  number  of  times  to  both  Q  and  P,  and  then 

match  the  resulting  cubic  polynomials.    Moreover,  once  the  relative  rotational  positions  of  P 

and  Q  have  been  estimated,  their  relative  translational  position  can  be  estimated  even  more 

easily  by  matching  the  terms  of  next-to-highest  degree  of  P  to  those  of  Q. 


-29- 

Each  of  the  coefficients  of  the  polynomial  Q  appears  quadratically  in  the  sum  (1')  to  be 
minimized.  If  the  rotated  polynomial  Q  is  known  to  belong  to  some  linear  subspace  of 
polynomials  not  necessarily  passing  through  the  origin  (or,  equivalently,  if  P  is  known  to 
belong  to  some  rotationally  invariant  linear  subspace  of  polynomials)  then  it  is  appropriate  to 
minimize  over  this  space  only.  Similarly  if  P  is  known  to  belong  to  some  rotationally 
invariant  quadratic  surface  in  the  space  of  polynomials,  it  is  appropriate  to  minimize  for  Q 
satisfying  this  constraint  also.  Beyond  this,  imposition  of  more  complex  constraints,  or  of 
more  than  one  quadratic  constraint,  significantly  increases  the  cost  of  minimization. 

If  the  polynomial  P  is  quadratic  it  will  satisfy  a  condition  AP  =  c  for  some  constant  c, 
and  since  this  linear  condition  is  rotationally  invariant  we  can  subject  Q  to  it  also.   Moreover, 

since  in  this  case  Pt= P  defines  a  vector  of  polynomials,  the  polynomial  B  =  '2lPl2  is  also 

dx, 

of  second  order  and  will  satisfy  a  second  condition  AS  =  c2,  which  is  quadratic  in  the 

coefficients  of  the  the  original  polynomial  P.  Hence  we  can  assume  that  Q  satisfies  this  same 

quadratic  condition.  It  follows  that  an  appropriate  way  of  matching  a  given  model  polynomial 

P  is  to  find  the  Q  which  minimizes 


dx, 


2 


{2Ax  +  v)(2Ar  +  v)  =  4x  A2x  +  4/Lc-  v  +  v2,     so    A^ 


dx, 


min  5>/  (C(x;))2,  (1») 

where  the  coefficients  c:  and  c2  are  taken  from  the  model  polynomial  P.  Note  that  if 
Q{x)  -  xAx  +  vx  +  b        is        a       quadratic       polynomial,        then       A£)  =  2tr  (A)        and 

=  8*r(A2). 

The  constraints  appearing  in  the  minimum  (1)  therefore  restrict  the  2a»  an<^  2a/y>  where 
A  =  \atj)  is  the  symmetric  matrix  of  the  quadratic  (i.e.  highest  order)  terms  in  the 
polynomial  Q. 

In  somewhat  more  detail:    The  sum  2w/(Q(*f)J     appearing  in  (1")  can  easily  be 
expressed  in  terms  of  the  coefficients  of  the  polynomial 
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Q(*)  =  '2,Cklmukvtz», 

where  x  =  («,  v,  z)  namely 

2  *i  (GW)2  =  2  era,  [s  *,  «f +  *  v,' +  <~  wr +  »j  i 

=  2  cWm  CI7,*  Kk(m,nm  I 

where  the  symmetric  matrix  AT  is  positive  definite.  (Note  that  each  of  the  entries  of  this 
matrix  is  simply  a  particular  weighted  moment  of  the  collection  of  surface  points  observed, 
this  moment  having  a  degree  which  is  at  most  twice  the  degree  of  the  polynomial  equation 
Q(x)  =  0  which  is  to  be  matched  to  these  surface  points.)  In  particular,  if 
Q(x)  =  x  A  x  +  b  x  +  c  is  a  quadratic  polynomial,  and  if  we  group  the  three  components 
of  b  together  with  the  scalar  c  to  form  a  4-dimensional  quantity  B  =  (b,  c),  we  can  write  the 
quadratic  form  (2)  as 

AKiA  +  2BK2A  +  BKyB,  (2*) 

and  so  can  rewrite  the  minimum  (1")  in  the  form 

min    ndn(AK1A  +  2BKJL  +  BK&)  (2") 

<r(A3)-c} 

=      min    A(K-K'K^K,)A. 

lr(A)-«r, 
1 

We  can  then  write  A  =  —cJ+A0,  where  tr(A0)  =  0,   allowing  the  minimum  (2")   to  be 


rewritten  as 


min         AKiA  +  ^-cAK.I+l-cllK.I-  (3) 

tr(A)-0  3     i  9 

fr(A2)-c,-cJ/3 


where 


i»^    ~    A|  /Li    ^3         Ai  . 

To  see  how  we  can  find  this  last  minimum,  let  the  variable  v  vary  over  the  space  of  all 
matrices  of  zero  trace,  and  let  v0  designate  the  matrix  c{  Af4  1/3.  Let  c3  =  c2  -  c\/3  and 
ct  =  cx2  IKt  1/9.  Then  the  minimum  (3)  can  be  written  as 
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min    vAf4v+2v0-v  +  c4  (4) 

By  Lagrange's  multiplier  principle  the  minimizing  vector  v  satisfies 
KAv  +  v0— \lv  =  0  for  some  real  X,  in  addition  to  the  equation  |v|2  =  cy  Thus  we 
have  v  =  (XZ-JSTJ-1^.  (However,  if  v0  is  orthogonal  to  the  eigenvector  c  of  £4  belonging  to 
a  particular  eigenvalue  u-  (in  particular,  if  v0  =  0,  which  will  always  be  the  case  if  c:  =  0)  we 
could  also  have  v  =  (p,-/f4)_1  va  +  ae  where  a  must  be  chosen  to  make  |v|2  =  c3.) 

Thus  (assuming  that  the  eigenvalues  of  the  positive  definite  symmetric  matrix  AT4  are 
simple)  we  need  only  examine  a  finite  number  of  possibilities  in  searching  for  the  v  which 
realizes  the  minimum  (4):  namely  these  cases  connected  with  eigenvalues  p.  of  KA  noted  just 
above,  plus  the  roots  X  of  the  equation  |(X7  -  ATJ-1  v0  \2  =  c3.  This  last  equation  is 
algebraic  and  of  degree  equal  to  the  dimension  of  the  space  3  by  3  of  matrices  of  trace  zero, 
i.e.  is  of  degree  10. 

Little  change  in  the  above  discussion  is  needed  to  cover  the  case  in  which  P  (hence  Q)  is 
of  even  degree  2*+ 2.   Here  we  have  only  to  replace  (1")  by 

min  2>,(q(x,))2  (1'") 

Introducing  Qk  =  A*  Q,  so  that  Ql  is  simply  a  quadratic  polynomial,  we  can  write  this  last 
minimum  as 


Afl, 


min  min     S^/fcW)2.  (5) 

The  'inner'  minimum  in  (5)  has  the  abstract  form 

min      qKq,  (6) 

Tq—a 

where  T  is  a  linear  transformation  A*  from  the  space  of  polynomials  of  order  2£  +  2  to  the 
space  of  polynomials  of  order  2  ,  and  AT  is  a  positive  definite  matrix.  Note  that  T  maps  its 
domain  onto  its  range;  to  see  this,  we  have  only  to  note  that  the  polynomial  x2,  whose  rotated 
versions  generate  the  whole  space  of  quadratic  polynomials,  belongs  to  the  range  of  A*,  as  is 
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clear  since  A*  x2k+2=  y  (2*  +  2)!  x2.  Putting  K^q  =  r  in  (6),  we  can  rewrite  (6)  as 

min     \r\2=   min  |r|2  (7) 

where  Tx  =  TK~1'2.  The  minimizing  r  in  (7)  is  clearly  the  unique  solution  of  Tx  r  =  a  which 
is  orthogonal  to  the  space  of  all  vectors  v  satisfying  7^  v  =  0,  We  can  express  this  v  as 
t  =  T\  {TxT\)~la,  since  then  we  have  Tx  r  =  (7^  7"*)  {TJ\)~l  a  =  a,  and  moreover  if 
7^  =  0  we  have  rv  =  T\  (Ji  Tj)'1  a-v  =  (Tx  T\Yl  a-Ti  v  =  0.  It  follows  that  the 
minimizing  q  in  (6)  is  AT1'2  r=  R-^  K'^2  T'1  CTK"1!0)"1  a  =  K~l  T\TK~l  T'Y^a, 
recovering  the  well  known  fact  that  the  minimum  (6)  is  just 

a  (TK-^T')'1  TK-lKK~lT"  (TK~l  T")-1  a  =  a  (TK~l  T)'1  a. 
This  shows  that  the  minimum  (5)  can  be  expressed  in  terms  of  the  polynomial  Q{  by 

minimum  of  the  form  that  we  have  already  treated,  namely 

min        Qi*iGi,  (8) 

AC,-*; 

where  the  symmetric  matrix  K  is  obtained  in  an  elementary  way  from  the  matrix  of  the 
quadratic  form  appearing  in  (5),  namely  by  two  matrix  inversions  and  two  matrix 
multiplications.   Finally,  the  minimum  (8)  can  be  found  in  the  manner  already  explained. 

When  we  apply  the  technique  just  described  to  determine  the  position  of  an  object  O 
several  of  whose  surfaces  may  be  visible,  it  is  essential  that  we  avoid  mixing  observed  points 
drawn  from  one  of  O's  surfaces  with  points  drawn  from  some  other  surface  of  O.  If  O's 
various  surfaces  are  separated  cleanly  enough  by  edges  at  which  the  tangent  surface  to  O 
changes  sharply,  it  may  be  possible  to  do  this  by  applying  an  edge  operator  to  an  intensity 
image  superimposed  on  a  depth  image,  or  directly  to  a  depth  image  to  bring  out  'jump' 
boundaries.  (Only  components  not  broken  by  edges  of  either  of  those  types  should  be  used 
for  matching.  Note  also  that  the  matching  method  that  we  have  described  is  likely  to  give 
immediate  warning  of  attempts  to  match  such  'mixed  observations',  since  such  attempts  are 
not  likely  to  generate  small  minima  (5),  given  that  the  two  parameters  cl  c2  appearing  in  (5) 


-33- 

are  both  taken  directly  from  one  of  the  finite  collection  of  model  surfaces  which  would  match 
the  scene  perfectly).  However,  cases  in  which  there  are  smooth  transitions  between  surfaces 
of  O  described  by  spline  surfaces,  satisfying  very  different  equations  but  constructed  to 
achieve  smooth  joins,  will  of  course  not  be  separable  by  any  visible  edges.  In  such  cases  the 
following  approach  may  be  useful.  Let  one  surface  satisfy  a  polynomial  equation  Pi  (x)  =  0, 
while  the  second  surface  satisfies  P2  (x)  =  0.  Then  the  union  of  the  two  surfaces  satisfies 
P:  P2  (x)  =  0,  and  can  therefore  be  matched  to  the  higher-degree  polynomial  Pl  P2.  A 
similar  observation  applies  to  corners  at  which  three  or  four  surface  patches  come  together. 

Sets  of  observed  points  for  use  in  the  matching  procedure  that  has  been  described  can 
be  gathered  in  any  convenient  way.  For  example,  one  can  gather  partial  depth  information 
by  illuminating  a  scene  with  just  one  plane  of  light,  which  gives  the  depth  of  all  object  points 
lying  in  the  illuminated  plane.  The  resulting  curve  can  broken  at  all  its  sharp  (particularly 
concave)  corners,  and  again  at  all  sharp  jumps  in  intensity,  yielding  curve  sections  each  of 
which  probably  lies  on  one  single  object  surface.  This  same  procedure  can  then  be  repeated 
for  a  second  plane  orthogonal  to  the  original  plane  of  illumination  and  cutting  it  at  a  point 
well  interior  to  one  of  the  aforementioned  curve  sections,  or,  still  better,  repeated  for  planes 
forming  a  pair  of  orthogonal  grids.  Any  two  3-D  curve  sections  likely  to  lie  together  in  one 
single  object  surface  then  'predict*  the  orientation  of  such  a  surface. 

Once  the  presence  and  orientation  of  a  particular  object  surface  has  been  predicted,  the 
position  in  space  of  the  body  containing  this  surface  can  be  determined,  and  then  the 
prediction  can  be  checked  by  comparing  an  artificially  generated  depth  image  of  the  body  to 
the  observed  scene  and  looking  for  substantial  regions  of  agreement.  Such  regions  should  of 
course  represent  other  surfaces  of  the  same  body.  (It  can  also  be  verified  that  no  point  in  the 
scene  actually  observed  lies  substantially  behind  any  surface  of  a  body  conjectured  to  be 
present).  Once  all  regions  of  an  observed  scene  have  been  identified  as  object  surfaces  in 
this  manner,  all  objects  present  and  visible  in  the  scene  will  be  known. 
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At  this  point,  the  further  task  of  following  the  motion  of  these  objects  becomes 
substantially  easier,  and  can  be  handled  as  follows:  for  each  visible  body,  one  can  take  points 
well  within  the  visible  surface  portions  P(x)  =  0.  Suppose  that  the  initial  position  of  each 
such  portion  P(x)  =  0  is  known,  and  suppose  that  after  a  small  interval  of  time  each  of  the 
points  of  this  surface  portion  has  moved  from  an  original  position  x  to  a  new  position  Rx+b, 
where  R=I  +  a  is  a  small  rotation  and  b  is  a  small  translation.   Suppose  that  points  x,  on  the 

body  surface  in  its  new  position  are  observed.  Then  the  points  /?-1  x,  —  b  lie  on  the  original 

i 
body  surface.    For  each  such  point,  take  any  sufficiently  close  point  y,  on  the  original  body 

surface  (as  known  from  its  equation  as  first  estimated  ),  for  example  let  v,  be  the  point  on 

the  original  surface  whose  projections  along  two  coordinate  axes  agree  with  those  of  x,;  these 

axes  must  be  chosen  so  that  the  third  coordinate  axis  intersects  the  surface  P(x)  =  0  at  a 

substantial  angle  near  x,.  Then  the  difference  vector  R~l  x,  —  b  —  y,  =  x,  —  y,  +  ax,  —  b  lies 

close  to  the  plane  tangent  to  P(x)  =  0  at  y,  ,  and  hence  can  be  taken  to  satisfy  the  equation 

Li  (•*/  -  V;  +  ax,  -  b)  =  0, 
where  L,  =  dP(y,)  is  the  vector  normal  to  P(x)  =  0  at  y,.  These  relationships  give  us  a  family 

of  equations  of  which  any  six  should  suffice  to  determine  the  six  linearly  independent 

quantities  [a ,b];  of  course,  it  is  better  to  use  more  equations  and  then  to  estimate  a  and  b 

using  a  least  squares  procedure.   Note  that  much  this  same  procedure  can  be  used  to  improve 

the  initially  conjectured  match  of  an  observed  body  to  a  model  body  once  the  observed 

surfaces    which    should   lie    in    each   particular   rotated/translated   model    body   have    been 

identified.    Note  also  that  if  the  manner  in  which  a  body's  shape  can  change  is  simple  and 

known  (as  in  the  case  of  an  expanding  sphere,  or  of  an  ellipsoid  gradually  changing  in 

eccentricity),   a  similar  simple  technique  can  be  used  to  track  shape  changes  as  well  as 

rotations. 

It  is  also  worth  observing  that  the  procedure  that  we  have  described  applies  individually 
to  each  of  the  observed  surfaces  in  a  scene  containing  multiple  bodies,  provided  only  that  the 
general  form  of  these  surfaces  is  known  (e.g.  it  may  be  known  that  every  such  surface  is 
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either  plane  or  quadric),  even  if  the  overall  way  in  which  these  surfaces  cohere  to  form 
bodies  is  not  known  initially.  Suppose,  for  example,  that  the  observed  scene  is  known  to 
contain  only  one  body.  Then  by  matching  the  surfaces  visible  from  one  viewpoint,  we  can 
define  part  of  the  body  shape,  and  by  doing  this  from  sufficiently  many  viewpoints  can  hope 
to  acquire  a  complete  understanding  of  the  object.  The  resulting  object  description  would 
consist  of  a  list  of  equations  P,  (x)  =  0  (or  more  precisely  of  inequalities  Pt  (x)a:  0)  defining 
all  the  surfaces  of  the  body  (the  side  of  these  surfaces  on  which  the  interior  of  the  body  lies 
being  described  by  the  aforementioned  inequalities),  together  with  a  set  of  lists,  each  of 
which  defines  the  body  surfaces  which  come  together  at  one  of  the  vertices  v  of  the  body,  in 
the  circular  order  of  these  surfaces  around  the  vertex  v. 

Development  and  demonstration  of  a  capability  to  acquire  the  shapes  of  manufactured 
objects  directly  from  combined  depth  and  intensity  image  observations  would  represent  a 
significant  step  forward  in  computer  vision  technology.  ° 
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